@---------------------------------------------@
@ GAUSS code for panel model with IID errors  @
@---------------------------------------------@

new;
cls;
library optmum;

@ ++++++++++++++++++++++++++++++++ Data +++++++++++++++++++++++++++++++++++++ @

load microdata[31631,4]=Micro.txt;
k=rows(microdata); 

/* Calculate growth rates and number of cross sections */

count=1;
datamid=zeros(1,5);

i=1;
do until i>(k-1);

z1=microdata[i+1,1];
z2=microdata[i+1,2];
z3=microdata[i+1,3]-microdata[i,3];
z4=microdata[i,3];
z5=microdata[i,4];
z=z1~z2~z3~z4~z5;

if microdata[i+1,1]==microdata[i,1];datamid=datamid|z;endif;
if microdata[i+1,1]>microdata[i,1];count=count+1;endif;

i=i+1;
endo;

kk=rows(datamid);
datamid=datamid[2:kk,.];
kkk=rows(datamid);

/* Calculate number of periods for every cross section */

t=ones(2069,1);
j=1;

i=1; do until i>(kkk-1);

if datamid[i+1,1]>datamid[i,1]; j=j+1;endif;
if datamid[i+1,1]==datamid[i,1]; t[j,1]=t[j,1]+1;endif;

i=i+1;
endo;

/* Datamid sorted by time */

datamid2=sortc(datamid,2);

/* Calculate number of cross sections for every period */

n=ones(25,1);
j=1;

i=1; do until i>(kkk-1);

if datamid2[i+1,2]>datamid2[i,2]; j=j+1;endif;
if datamid2[i+1,2]==datamid2[i,2]; n[j,1]=n[j,1]+1;endif;

i=i+1;
endo;

/* Calculate cross section fixed effects */

lambdai=zeros(2069,1);
sto=0;

j=1;
do until j>2069;

sta=sto+1;
sto=sto+t[j,1];

lambdai[j,1]=meanc(datamid[sta:sto,3]);

j=j+1;
endo;

/* Calculate time fixed effects */

lambdat=zeros(25,1);
sto=0;

j=1;
do until j>25;

sta=sto+1;
sto=sto+n[j,1];

lambdat[j,1]=meanc(datamid2[sta:sto,3]);

j=j+1;
endo;

/* Calculate transformed growth rates */

datafinal=zeros(kkk,4);
datafinal[.,1:2]=datamid[.,1:2];
year=datamid[1,2];
datafinal[1,3]=datamid[1,3]-lambdat[year-67,1];
datafinal[1,4]=lambdat[year-67,1];

i=1;do until i>(kkk-1);

year=datamid[i+1,2];
datafinal[i+1,3]=datamid[i+1,3]-lambdat[year-67,1];
datafinal[i+1,4]=lambdat[year-67,1];

i=i+1;
endo;

@ ++++++++++++++++++++++++++++++++ Estimation +++++++++++++++++++++++++++++++++++++ @

par=2;     
par_in1=-2|-1;
{xout,fout,gout,cout}= optmum(&li,par_in1);
par_final=trans(xout);
teta=par_final[1,1];  
sigma=par_final[2,1];

if cout==0;
LL_val=-fout;
hessian = hessp(&li,xout);
if minc(eigrs(hessian))<=0.00000001;
cls;
print "   ";
print "Numerical Hessian using 'hessp' is not positive definite, using '_opfhess' instead";
omega_hat=inv(_opfhess);
else;
cls;
omega_hat=inv(hessian);
endif;
endif;

grdn=gradp(&trans,xout); 
cov=grdn*omega_hat*grdn';     
sd =((sqrt(diag(cov))));
AIC=(-2*LL_val+2*par)/kkk;

@ +++++++++++++++++++++++++++++++ Output ++++++++++++++++++++++++++++++++++++++ @

gmm=zeros(kkk,8);
gmm[.,1:2]=datafinal[.,1:2];
gmm[.,3]=datamid[.,5];
gmm[.,4]=datamid[.,4];

gmm=fin(par_final);
gmm2=sortc(gmm,2);

var=zeros(25,1);

sto=0;
j=1;
do until j>25;

sta=sto+1;
sto=sto+n[j,1];

var[j,1]=meanc(gmm2[sta:sto,6]);

j=j+1;
endo;

sqnphi=sqrt(n)~sqrt(n);
par_phi=teta|sigma;
alfa=0;beta=0;gama=0;

@ ++++++++++++++++++++++++++++++++ Report +++++++++++++++++++++++++++++++++++++ @

print "  ";
print "********************************** Estimation Output *****************************************";
print " ";
print "MA(1) parameter and unconditional variance shock (left) and standard errors (right)";
print par_final~sd ;
print "**********************************************************************************************";
print "  ";
print "Maximized log-likelihood (left) and Akaike information criterion (right)";
print LL_val~AIC;
print "  ";

@ ++++++++++++++++++++++++++++++++ Likelihood +++++++++++++++++++++++++++++++++++++ @

proc li(prmtr);
local par_in,teta,vareps,sum,sumli,sto,sta;
local i,j,Y,vt,vt2,vart;
 
par_in=trans(prmtr);
teta=par_in[1,1];
vareps=par_in[2,1];

sum=0;sto=0;

j=1;
do until j>count;

sta=sto+1;
sto=sto+t[j,1];

Y=datafinal[sta:sto,3];
vt=0;
vt2=vareps;

i=1; do until i>t[j,1];

vart=vareps;
vt=Y[i,1]+teta*vt;
vt2=vt^2;

sum=sum+ln(vart)+vt2*inv(vart);

i=i+1; endo;

j=j+1;
endo;

sumli=-(-0.5*(kkk*ln(2*pi))-0.5*sum);

retp(sumli);
endp;

@ ++++++++++++++++++++++++++++++++ Data +++++++++++++++++++++++++++++++++++++ @

proc fin(prmtr);
local par_in,teta,vareps,sto,sta;
local i,j,Y,vt,vt2,vart,shock,variance;
 
par_in=prmtr;
teta=par_in[1,1];
vareps=par_in[2,1];

sto=0;

j=1;
do until j>count;

sta=sto+1;
sto=sto+t[j,1];

Y=datafinal[sta:sto,3];
vt=0;
vt2=vareps;
shock=zeros(t[j,1],1);
variance=zeros(t[j,1],1);

i=1; do until i>t[j,1];

vart=vareps;
variance[i,1]=vart;
vt=Y[i,1]+teta*vt;
shock[i,1]=vt;
vt2=vt^2;

i=i+1; endo;

gmm[sta:sto,5]=shock;
gmm[sta:sto,6]=variance;

j=j+1;
endo;

retp(gmm);
endp;

@ ++++++++++++++++++++++++++++++++ Transformation +++++++++++++++++++++++++++++++++++++ @

proc trans(x0); 
local x1;

x1=x0;

x1[1,1]=1/(1+exp(-x0[1,1]));  
x1[2,1]=exp(x0[2,1]);

retp(x1);
endp;

